Fungal OTU’s across life stages of boreal toad: published and supplemental analyses
Alexandra Alexiev
12/17/2019
#packages and set working directory
library(vegan)
library(ggplot2)
library(mctoolsr)
library(tidyr)
library(tibble)
library(gridExtra)
library(dplyr)
library(ggpubr)
library(cowplot)
library(kableExtra)
library(RColorBrewer)
set.seed(33333)
work_dir <- "/Users/alal3493/Documents/Projects/DevelopmentalBorealToad/FungalDevData/02_PublishedAnalysis/"
seq_dir <- "/Users/alal3493/Documents/Projects/DevelopmentalBorealToad/FungalDevData/00_SeqProcessing/"Basic measures of highly abundant or present taxa
#load in starting files
mapfp_all <- paste0(work_dir,"Allsamps_clean_maptab.txt")
tabfp_all <- paste0(work_dir,"Allsamps_clean_otutab.txt")
input_all <- load_taxa_table(tab_fp = tabfp_all, map_fp = mapfp_all) # this file has both toad and environment sample types in it (172 samples)
mapfp_t <- paste0(work_dir,"Toad_clean_maptab.txt")
tabfp_t <- paste0(work_dir,"Toad_clean_otutab.txt")
input_toad <- load_taxa_table(tab_fp = tabfp_t, map_fp = mapfp_t) # this file only has toad samples in it (124 samples)
# order by life stage
life_list <- c("Soil", "Sediment", "Water",
"Eggs.11.15", "Tadpole.20.22",
"Tadpole.23.25", "Tadpole.25.27",
"Tadpole.27.29", "Tadpole.29.31",
"Tadpole.31.35", "Tadpole.36.39",
"Metamorph.40.46", "Subadult", "Adult")
input_all$map_loaded$Gosner_Stage <-
factor(input_all$map_loaded$Gosner_Stage,
levels=life_list)
input_toad$map_loaded$Gosner_Stage <-
factor(input_toad$map_loaded$Gosner_Stage,
levels=life_list)
# and order the other column I'll be using as well, where tadpoles are lumped
samptype_list <- c("Soil", "Sediment", "Water", "Eggs",
"Tadpole", "Metamorph", "Subadult", "Adult")
input_all$map_loaded$Life_Stage_Simplified <-
factor(input_all$map_loaded$Life_Stage_Simplified,
levels=samptype_list)
input_toad$map_loaded$Life_Stage_Simplified <-
factor(input_toad$map_loaded$Life_Stage_Simplified,
levels=samptype_list)# number of total fungal taxa found in all samples
nrow(input_all$data_loaded) # 2243
## [1] 2243
nrow(input_toad$data_loaded) # 1703
## [1] 1760
# top ten taxa on toads and in environment
toptaxa_all <- return_top_taxa(input = input_all, number_taxa = 10)
toptaxa_toad <- return_top_taxa(input = input_toad, number_taxa = 10)
kable(toptaxa_all, caption="Top ten taxa in toad and environment samples") %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| taxonomy1 | taxonomy2 | taxonomy3 | taxonomy4 | taxonomy5 | taxonomy6 | taxonomy7 | |
|---|---|---|---|---|---|---|---|
| OTU_1 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Capnodiales | f__Mycosphaerellaceae | g__Mycosphaerella | s__Mycosphaerella_tassiana |
| OTU_2 | k__Fungi | NA | NA | NA | NA | NA | NA |
| OTU_4 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Capnodiales | f__Cladosporiaceae | g__Cladosporium | s__Cladosporium_delicatulum |
| OTU_3 | k__Fungi | p__Ascomycota | c__Leotiomycetes | o__Helotiales | f__Hyaloscyphaceae | g__Pezizella | s__Pezizella_epithallina |
| OTU_7180 | k__Fungi | p__Basidiomycota | c__Agaricomycetes | o__Agaricales | f__Schizophyllaceae | g__Schizophyllum | s__Schizophyllum_commune |
| OTU_34 | k__Fungi | p__Rozellomycota | c__unidentified | o__unidentified | f__unidentified | g__unidentified | s__unidentified |
| OTU_8 | k__Fungi | p__Ascomycota | c__Leotiomycetes | o__Thelebolales | f__Pseudeurotiaceae | g__Pseudeurotium | s__Pseudeurotium_hygrophilum |
| OTU_11 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Pleosporales | f__Phaeosphaeriaceae | g__unidentified | s__unidentified |
| OTU_37 | k__Fungi | p__Basidiomycota | c__Tremellomycetes | o__Filobasidiales | f__Filobasidiaceae | g__Naganishia | s__Naganishia_albida |
| OTU_15 | k__Fungi | p__Basidiomycota | c__Tremellomycetes | o__Filobasidiales | f__unidentified | g__unidentified | s__unidentified |
kable(toptaxa_toad, caption="Top ten taxa in toad samples") %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| taxonomy1 | taxonomy2 | taxonomy3 | taxonomy4 | taxonomy5 | taxonomy6 | taxonomy7 | |
|---|---|---|---|---|---|---|---|
| OTU_1 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Capnodiales | f__Mycosphaerellaceae | g__Mycosphaerella | s__Mycosphaerella_tassiana |
| OTU_4 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Capnodiales | f__Cladosporiaceae | g__Cladosporium | s__Cladosporium_delicatulum |
| OTU_2 | k__Fungi | NA | NA | NA | NA | NA | NA |
| OTU_7180 | k__Fungi | p__Basidiomycota | c__Agaricomycetes | o__Agaricales | f__Schizophyllaceae | g__Schizophyllum | s__Schizophyllum_commune |
| OTU_11 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Pleosporales | f__Phaeosphaeriaceae | g__unidentified | s__unidentified |
| OTU_3 | k__Fungi | p__Ascomycota | c__Leotiomycetes | o__Helotiales | f__Hyaloscyphaceae | g__Pezizella | s__Pezizella_epithallina |
| OTU_15 | k__Fungi | p__Basidiomycota | c__Tremellomycetes | o__Filobasidiales | f__unidentified | g__unidentified | s__unidentified |
| OTU_34 | k__Fungi | p__Rozellomycota | c__unidentified | o__unidentified | f__unidentified | g__unidentified | s__unidentified |
| OTU_37 | k__Fungi | p__Basidiomycota | c__Tremellomycetes | o__Filobasidiales | f__Filobasidiaceae | g__Naganishia | s__Naganishia_albida |
| OTU_12 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Pleosporales | f__Phaeosphaeriaceae | g__Sclerostagonospora | s__Sclerostagonospora_phragmiticola |
Q1: What are the factors associated with variation in fungal communities?
Alpha diversity differences between life stages and environment samples
Below is a supplementary figure comparing three alpha diversity measures across all the sample types. Since site is a confounding factor here, given exploratory analyses, I analyzed with blocking by site (“strata”).
# calculate alpha diversity measures
input_all$map_loaded$rich = calc_diversity(
tax_table = input_all$data_loaded, metric = 'richness')
input_all$map_loaded$simpson = calc_diversity(
tax_table = input_all$data_loaded, metric = 'simpson')
input_all$map_loaded$shannon = calc_diversity(
tax_table = input_all$data_loaded, metric = 'shannon')
# Kruskall-Wallis test on life stages, strata-site
# species richness
krusk_rich <- kruskal.test(input_all$map_loaded$rich ~
input_all$map_loaded$Gosner_Stage)
bonpval_rich <- p.adjust(p = krusk_rich$p.value, method = "bonferroni",
n = length(unique(input_all$map_loaded$Gosner_Stage)))
krusk_rich; bonpval_rich
##
## Kruskal-Wallis rank sum test
##
## data: input_all$map_loaded$rich by input_all$map_loaded$Gosner_Stage
## Kruskal-Wallis chi-squared = 106.55, df = 13, p-value < 2.2e-16
## [1] 1.235034e-15
# simpson (evenness)
krusk_simp <- kruskal.test(input_all$map_loaded$simpson ~
input_all$map_loaded$Gosner_Stage)
bonpval_simp <- p.adjust(p = krusk_simp$p.value, method = "bonferroni",
n = length(unique(input_all$map_loaded$Gosner_Stage)))
krusk_simp; bonpval_simp
##
## Kruskal-Wallis rank sum test
##
## data: input_all$map_loaded$simpson by input_all$map_loaded$Gosner_Stage
## Kruskal-Wallis chi-squared = 29.009, df = 13, p-value = 0.006528
## [1] 0.09138674
# shannon (both richness and evenness)
krusk_shan <- kruskal.test(input_all$map_loaded$shannon ~
input_all$map_loaded$Gosner_Stage)
bonpval_shan <- p.adjust(p = krusk_shan$p.value, method = "bonferroni",
n = length(unique(input_all$map_loaded$Gosner_Stage)))
krusk_shan; bonpval_shan
##
## Kruskal-Wallis rank sum test
##
## data: input_all$map_loaded$shannon by input_all$map_loaded$Gosner_Stage
## Kruskal-Wallis chi-squared = 52.554, df = 13, p-value = 1.08e-06
## [1] 1.51244e-05
# GRAPHING
# Species richness graphs: mean number of OTUs found per sample type
liferich <- ggplot(input_all$map_loaded, aes(Gosner_Stage, rich)) +
geom_boxplot() + xlab("") +
ylab("Fungal OTU Richness") + ggtitle("A") +
geom_text(x=7,y=194, size=3,
label=paste0("KW chi-sq stat=",signif(krusk_rich$statistic, digits=5),
"\n","p=",signif(bonpval_rich, digits=3))) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
# Simpson's: measures species/OTU evenness, or how close in number each species in the environment
lifesimp <- ggplot(input_all$map_loaded, aes(Gosner_Stage, simpson)) +
geom_boxplot() + xlab("Sample Type") +
ylab("Simpson's Diversity Metric") + ylim(0,1.2) + ggtitle("B") +
geom_text(x=7,y=1.15, size=3,
label=paste0("KW chi-sq stat=",signif(krusk_simp$statistic, digits=5),
"\n","p=",signif(bonpval_simp, digits=3))) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
# Shannon: combination of richness and evenness
lifeshan <- ggplot(input_all$map_loaded, aes(Gosner_Stage, shannon)) +
geom_boxplot() + xlab("") +
ylab("Shannon Diversity Index") + ggtitle("C") +
geom_text(x=7,y=4.1, size=3,
label=paste0("KW chi-sq stat=",signif(krusk_shan$statistic, digits=5),
"\n","p=",signif(bonpval_shan, digits=3))) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
alpha_panel <- grid.arrange(liferich, lifesimp, lifeshan, ncol = 3)
alpha_panel
## TableGrob (1 x 3) "arrange": 3 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (1-1,3-3) arrange gtable[layout]
# ggsave(paste0(work_dir,"Figs_Tables/AllAlphaDiv.png"), alpha_panel, dpi=300,
# width=10, height=7)Pairwise OTU richness Stats
- Here, I’m using a two-sample Mann-Whitney test (or unpaired wilcoxon test). This is a non-parametric test of whether two populations are different across two different life stages, while still assuming those populations are independent. Since my samples came from different animals at different life stages, I need to make this assumption, so I cannot use a paired wilcoxon test.
- I calculated Bonferroni (signif if p<=0.05) and FDR (p<=0.05) p-value corrections, used both together to figure out how big of a difference there was between each life stage.
# Are pairwise alpha diversities same or different across life stages?
nlevels(input_all$map_loaded$Gosner_Stage) # 14
## [1] 14
cats <- levels(input_all$map_loaded$Gosner_Stage) # define category levels
cats_pairs <- combn(cats, 2) # make pairs of these
pvals <- c() # output df for pvalues
Wstat <- c() # output df for Wstat
for ( i in 1:ncol(cats_pairs) ) { # for each pair of sample types
pair <- cats_pairs[, i] # define the pair
Gosner_rich <- input_all$map_loaded %>% # take the mapping file
dplyr::select(Gosner_Stage, rich) %>% # select these two columns from it
filter(Gosner_Stage %in% pair) # filter the sample types in the pair
wil <- wilcox.test(rich ~ Gosner_Stage, data = Gosner_rich, paired = F, exact = F) # run Wilcoxon test
pvals <- c(pvals, wil$p.value) # save the p val
Wstat <- c(Wstat, wil$statistic) # save the W statistic
}
results <- data.frame(t(cats_pairs), Wstat, pvals) # make a dataframe of the useful outputs
results$pvalBon = pvals * length(pvals) # calc bonferroni pval correction
results$pvalFDR = round(pvals * (length(pvals) / rank(pvals, ties.method = "average")),
3) # calc FDR pval correction
results <- results %>%
arrange(pvalBon, pvalFDR) %>% # arrange in ascending order by the corrected pvals
filter(pvalBon <= 0.05) %>%
filter(pvalFDR <= 0.05) %>% # filter by the two corrected pvals being <= 0.05
dplyr::select(-pvals) # remove pval column
# write.table(results, paste0(work_dir, "Figs_Tables/pairedMW_alphadiv.txt"),
# sep="\t", quote = F, row.names = F)Make three paneled figure with alpha and beta diversity measures
# make color palette
nlevels(input_all$map_loaded$Life_Stage_Simplified)
## [1] 8
brewer.pal(n = 8, name = "Dark2") # display hexcodes for the pallete I'm using but want to change slightly
## [1] "#1B9E77" "#D95F02" "#7570B3" "#E7298A" "#66A61E" "#E6AB02" "#A6761D"
## [8] "#666666"
samptype_cols <- c("#666666", "#A6761D",
"#3a6b94", "#66A61E",
"#E7298A", "#7570B3",
"#D95F02", "#1B9E77")
# make OTU richness figure
OTUrich <- ggplot(input_all$map_loaded, aes(Gosner_Stage, rich)) +
geom_boxplot(alpha=0.6, outlier.shape=NA, aes(fill=Life_Stage_Simplified)) + xlab("") +
ylab("Fungal OTU Richness") + ggtitle("A") +
geom_point(aes(fill=Life_Stage_Simplified), size=3, shape=21, position=position_jitter(
width=0.25
)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
axis.title = element_text(size=16),
axis.text = element_text(size=14),
legend.title = element_text(size=16),
legend.text = element_text(size=14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(size = 30),
legend.position="none") +
labs(fill="Sample Type") +
scale_fill_manual(values=samptype_cols) +
scale_x_discrete(labels=c("Soil", "Sediment", "Water",
"Eggs", "Tadpole (20-22)",
"Tadpole (23-25)", "Tadpole (25-27)",
"Tadpole (27-29)", "Tadpole (29-31)",
"Tadpole (31-35)", "Tadpole (36-39)",
"Metamorph (40-46)", "Subadult", "Adult"))
# make NMDS plot
# create distance matrix
md_transformed <- t(sqrt(input_all$data_loaded))
dm <- vegdist(md_transformed, method = "bray")
# run NMDS
md.nmds <- metaMDS(dm, k = 3, trymax = 1000) #solution reached, stress=0.19
## Run 0 stress 0.1991807
## Run 1 stress 0.198966
## ... New best solution
## ... Procrustes: rmse 0.03891056 max resid 0.2153907
## Run 2 stress 0.2031021
## Run 3 stress 0.19948
## Run 4 stress 0.2002888
## Run 5 stress 0.2003817
## Run 6 stress 0.1986503
## ... New best solution
## ... Procrustes: rmse 0.03456291 max resid 0.1671252
## Run 7 stress 0.1988183
## ... Procrustes: rmse 0.03830062 max resid 0.1667616
## Run 8 stress 0.1993512
## Run 9 stress 0.2004438
## Run 10 stress 0.1986031
## ... New best solution
## ... Procrustes: rmse 0.01465033 max resid 0.1623897
## Run 11 stress 0.1989448
## ... Procrustes: rmse 0.04406373 max resid 0.2326543
## Run 12 stress 0.1991588
## Run 13 stress 0.1996618
## Run 14 stress 0.2092705
## Run 15 stress 0.1989837
## ... Procrustes: rmse 0.04302669 max resid 0.2298482
## Run 16 stress 0.2000927
## Run 17 stress 0.1989927
## ... Procrustes: rmse 0.01564361 max resid 0.1366187
## Run 18 stress 0.2017775
## Run 19 stress 0.1992631
## Run 20 stress 0.2013856
## Run 21 stress 0.1998188
## Run 22 stress 0.1987841
## ... Procrustes: rmse 0.03513183 max resid 0.1730398
## Run 23 stress 0.198889
## ... Procrustes: rmse 0.03759162 max resid 0.2387053
## Run 24 stress 0.2008108
## Run 25 stress 0.1997658
## Run 26 stress 0.1987944
## ... Procrustes: rmse 0.03852934 max resid 0.1676246
## Run 27 stress 0.1999514
## Run 28 stress 0.1991686
## Run 29 stress 0.19848
## ... New best solution
## ... Procrustes: rmse 0.03077958 max resid 0.1809424
## Run 30 stress 0.1987633
## ... Procrustes: rmse 0.0287966 max resid 0.1546369
## Run 31 stress 0.1987755
## ... Procrustes: rmse 0.01226734 max resid 0.10798
## Run 32 stress 0.1990778
## Run 33 stress 0.1992444
## Run 34 stress 0.1997452
## Run 35 stress 0.1988952
## ... Procrustes: rmse 0.02977232 max resid 0.1597155
## Run 36 stress 0.1995298
## Run 37 stress 0.1996322
## Run 38 stress 0.1994623
## Run 39 stress 0.1987867
## ... Procrustes: rmse 0.0239836 max resid 0.2207002
## Run 40 stress 0.2007746
## Run 41 stress 0.1990494
## Run 42 stress 0.2000423
## Run 43 stress 0.2008447
## Run 44 stress 0.1992178
## Run 45 stress 0.2021716
## Run 46 stress 0.2001571
## Run 47 stress 0.1986133
## ... Procrustes: rmse 0.01528265 max resid 0.1140558
## Run 48 stress 0.1998996
## Run 49 stress 0.2015246
## Run 50 stress 0.1995329
## Run 51 stress 0.1995766
## Run 52 stress 0.1987726
## ... Procrustes: rmse 0.02923879 max resid 0.156082
## Run 53 stress 0.1988128
## ... Procrustes: rmse 0.01379297 max resid 0.1567985
## Run 54 stress 0.1987881
## ... Procrustes: rmse 0.01621759 max resid 0.1491172
## Run 55 stress 0.1995058
## Run 56 stress 0.1988069
## ... Procrustes: rmse 0.02723704 max resid 0.1558522
## Run 57 stress 0.1984979
## ... Procrustes: rmse 0.002864535 max resid 0.03240811
## Run 58 stress 0.1982229
## ... New best solution
## ... Procrustes: rmse 0.01053161 max resid 0.1069317
## Run 59 stress 0.2042082
## Run 60 stress 0.1988679
## Run 61 stress 0.1988407
## Run 62 stress 0.1999197
## Run 63 stress 0.2023767
## Run 64 stress 0.1986061
## ... Procrustes: rmse 0.02027476 max resid 0.15781
## Run 65 stress 0.2004686
## Run 66 stress 0.1986067
## ... Procrustes: rmse 0.03168572 max resid 0.1844575
## Run 67 stress 0.1993562
## Run 68 stress 0.1989221
## Run 69 stress 0.199448
## Run 70 stress 0.1987401
## Run 71 stress 0.1988833
## Run 72 stress 0.1986434
## ... Procrustes: rmse 0.03020556 max resid 0.1832149
## Run 73 stress 0.1993182
## Run 74 stress 0.1996823
## Run 75 stress 0.2019042
## Run 76 stress 0.1989363
## Run 77 stress 0.199273
## Run 78 stress 0.1987613
## Run 79 stress 0.1989813
## Run 80 stress 0.1990984
## Run 81 stress 0.201643
## Run 82 stress 0.2011198
## Run 83 stress 0.1988251
## Run 84 stress 0.1990841
## Run 85 stress 0.199456
## Run 86 stress 0.2003537
## Run 87 stress 0.2004725
## Run 88 stress 0.1996851
## Run 89 stress 0.1989552
## Run 90 stress 0.1992805
## Run 91 stress 0.1984808
## ... Procrustes: rmse 0.01049447 max resid 0.1069944
## Run 92 stress 0.1995911
## Run 93 stress 0.2009097
## Run 94 stress 0.2020523
## Run 95 stress 0.1989221
## Run 96 stress 0.1985056
## ... Procrustes: rmse 0.0308026 max resid 0.184883
## Run 97 stress 0.2006718
## Run 98 stress 0.1999769
## Run 99 stress 0.1997734
## Run 100 stress 0.1996066
## Run 101 stress 0.2008897
## Run 102 stress 0.1987969
## Run 103 stress 0.1997936
## Run 104 stress 0.1995283
## Run 105 stress 0.19823
## ... Procrustes: rmse 0.004058423 max resid 0.03202131
## Run 106 stress 0.198935
## Run 107 stress 0.1988275
## Run 108 stress 0.1988211
## Run 109 stress 0.2132471
## Run 110 stress 0.1996344
## Run 111 stress 0.1988408
## Run 112 stress 0.1989723
## Run 113 stress 0.1990989
## Run 114 stress 0.2014891
## Run 115 stress 0.2004305
## Run 116 stress 0.2013293
## Run 117 stress 0.1994356
## Run 118 stress 0.1996882
## Run 119 stress 0.20103
## Run 120 stress 0.1996745
## Run 121 stress 0.1989367
## Run 122 stress 0.1989288
## Run 123 stress 0.1987036
## ... Procrustes: rmse 0.03038804 max resid 0.1899616
## Run 124 stress 0.1991161
## Run 125 stress 0.1999662
## Run 126 stress 0.1986356
## ... Procrustes: rmse 0.01922748 max resid 0.2101567
## Run 127 stress 0.1995036
## Run 128 stress 0.2027145
## Run 129 stress 0.1993594
## Run 130 stress 0.2120321
## Run 131 stress 0.1997443
## Run 132 stress 0.1991384
## Run 133 stress 0.2018665
## Run 134 stress 0.1988933
## Run 135 stress 0.1985443
## ... Procrustes: rmse 0.01156856 max resid 0.0983231
## Run 136 stress 0.2002259
## Run 137 stress 0.199142
## Run 138 stress 0.1988366
## Run 139 stress 0.2047048
## Run 140 stress 0.1987794
## Run 141 stress 0.2019196
## Run 142 stress 0.1998272
## Run 143 stress 0.1997401
## Run 144 stress 0.2032514
## Run 145 stress 0.1996975
## Run 146 stress 0.1988886
## Run 147 stress 0.1988603
## Run 148 stress 0.1990889
## Run 149 stress 0.1997159
## Run 150 stress 0.1990913
## Run 151 stress 0.2001639
## Run 152 stress 0.1983768
## ... Procrustes: rmse 0.007243963 max resid 0.07257952
## Run 153 stress 0.1985728
## ... Procrustes: rmse 0.03036131 max resid 0.1843549
## Run 154 stress 0.1988148
## Run 155 stress 0.198794
## Run 156 stress 0.1989918
## Run 157 stress 0.1989525
## Run 158 stress 0.2087058
## Run 159 stress 0.2006631
## Run 160 stress 0.1990916
## Run 161 stress 0.199112
## Run 162 stress 0.1989392
## Run 163 stress 0.1995455
## Run 164 stress 0.1985628
## ... Procrustes: rmse 0.02544981 max resid 0.1691289
## Run 165 stress 0.1994956
## Run 166 stress 0.198505
## ... Procrustes: rmse 0.0150738 max resid 0.1561891
## Run 167 stress 0.1991978
## Run 168 stress 0.1993111
## Run 169 stress 0.1992908
## Run 170 stress 0.1984406
## ... Procrustes: rmse 0.01871372 max resid 0.1549253
## Run 171 stress 0.1987424
## Run 172 stress 0.198577
## ... Procrustes: rmse 0.03053336 max resid 0.1863732
## Run 173 stress 0.1987948
## Run 174 stress 0.2074785
## Run 175 stress 0.2032761
## Run 176 stress 0.1986036
## ... Procrustes: rmse 0.03159869 max resid 0.1840107
## Run 177 stress 0.1991972
## Run 178 stress 0.1989689
## Run 179 stress 0.2000027
## Run 180 stress 0.1988215
## Run 181 stress 0.1988701
## Run 182 stress 0.199022
## Run 183 stress 0.1995071
## Run 184 stress 0.2000039
## Run 185 stress 0.1982403
## ... Procrustes: rmse 0.003224357 max resid 0.03028165
## Run 186 stress 0.1984985
## ... Procrustes: rmse 0.006869526 max resid 0.08173671
## Run 187 stress 0.1988268
## Run 188 stress 0.1993313
## Run 189 stress 0.1993146
## Run 190 stress 0.1981851
## ... New best solution
## ... Procrustes: rmse 0.006471767 max resid 0.07292009
## Run 191 stress 0.1986486
## ... Procrustes: rmse 0.02142667 max resid 0.1580042
## Run 192 stress 0.1988318
## Run 193 stress 0.1989523
## Run 194 stress 0.1981788
## ... New best solution
## ... Procrustes: rmse 0.009680974 max resid 0.119532
## Run 195 stress 0.1994719
## Run 196 stress 0.2006912
## Run 197 stress 0.1994413
## Run 198 stress 0.2028805
## Run 199 stress 0.1990498
## Run 200 stress 0.2014256
## Run 201 stress 0.1995668
## Run 202 stress 0.1993928
## Run 203 stress 0.1989151
## Run 204 stress 0.1995219
## Run 205 stress 0.1985167
## ... Procrustes: rmse 0.02839478 max resid 0.1746835
## Run 206 stress 0.1987447
## Run 207 stress 0.1991223
## Run 208 stress 0.1992759
## Run 209 stress 0.1994279
## Run 210 stress 0.1988712
## Run 211 stress 0.1995471
## Run 212 stress 0.201507
## Run 213 stress 0.1990177
## Run 214 stress 0.2019303
## Run 215 stress 0.1987481
## Run 216 stress 0.2043019
## Run 217 stress 0.2117713
## Run 218 stress 0.2001185
## Run 219 stress 0.2017796
## Run 220 stress 0.2101078
## Run 221 stress 0.199161
## Run 222 stress 0.1991476
## Run 223 stress 0.2098873
## Run 224 stress 0.2002368
## Run 225 stress 0.1987793
## Run 226 stress 0.1997222
## Run 227 stress 0.199483
## Run 228 stress 0.1997935
## Run 229 stress 0.2000397
## Run 230 stress 0.199509
## Run 231 stress 0.2089883
## Run 232 stress 0.2012403
## Run 233 stress 0.199842
## Run 234 stress 0.199508
## Run 235 stress 0.2023137
## Run 236 stress 0.199538
## Run 237 stress 0.1991486
## Run 238 stress 0.2000331
## Run 239 stress 0.1986312
## ... Procrustes: rmse 0.01557071 max resid 0.1623782
## Run 240 stress 0.2045712
## Run 241 stress 0.2015893
## Run 242 stress 0.2032
## Run 243 stress 0.1999276
## Run 244 stress 0.1985062
## ... Procrustes: rmse 0.02855012 max resid 0.1765339
## Run 245 stress 0.1985494
## ... Procrustes: rmse 0.03131395 max resid 0.1803953
## Run 246 stress 0.1992457
## Run 247 stress 0.201894
## Run 248 stress 0.202881
## Run 249 stress 0.1987257
## Run 250 stress 0.2009174
## Run 251 stress 0.1990257
## Run 252 stress 0.1990017
## Run 253 stress 0.2139191
## Run 254 stress 0.1999195
## Run 255 stress 0.2014226
## Run 256 stress 0.1999163
## Run 257 stress 0.1990871
## Run 258 stress 0.1993242
## Run 259 stress 0.199032
## Run 260 stress 0.1995299
## Run 261 stress 0.2041143
## Run 262 stress 0.1993912
## Run 263 stress 0.1987431
## Run 264 stress 0.1988431
## Run 265 stress 0.1987841
## Run 266 stress 0.198789
## Run 267 stress 0.1986663
## ... Procrustes: rmse 0.02858086 max resid 0.1729659
## Run 268 stress 0.1983763
## ... Procrustes: rmse 0.01417369 max resid 0.1575988
## Run 269 stress 0.1994074
## Run 270 stress 0.199073
## Run 271 stress 0.2118861
## Run 272 stress 0.199112
## Run 273 stress 0.199603
## Run 274 stress 0.1988698
## Run 275 stress 0.1985094
## ... Procrustes: rmse 0.01214482 max resid 0.1115698
## Run 276 stress 0.1993897
## Run 277 stress 0.1985726
## ... Procrustes: rmse 0.02770376 max resid 0.1752979
## Run 278 stress 0.1994618
## Run 279 stress 0.199095
## Run 280 stress 0.1994405
## Run 281 stress 0.1991439
## Run 282 stress 0.1989474
## Run 283 stress 0.1985908
## ... Procrustes: rmse 0.0277561 max resid 0.1738694
## Run 284 stress 0.1995117
## Run 285 stress 0.1985605
## ... Procrustes: rmse 0.02361725 max resid 0.1675909
## Run 286 stress 0.1990911
## Run 287 stress 0.1994633
## Run 288 stress 0.201189
## Run 289 stress 0.203217
## Run 290 stress 0.2015123
## Run 291 stress 0.1986045
## ... Procrustes: rmse 0.02780458 max resid 0.1740357
## Run 292 stress 0.1996768
## Run 293 stress 0.1995264
## Run 294 stress 0.2024634
## Run 295 stress 0.1995038
## Run 296 stress 0.1991955
## Run 297 stress 0.1987386
## Run 298 stress 0.1992106
## Run 299 stress 0.1993784
## Run 300 stress 0.1994614
## Run 301 stress 0.2024684
## Run 302 stress 0.1995412
## Run 303 stress 0.2017279
## Run 304 stress 0.2015048
## Run 305 stress 0.1991423
## Run 306 stress 0.1995625
## Run 307 stress 0.1988322
## Run 308 stress 0.1991954
## Run 309 stress 0.1985902
## ... Procrustes: rmse 0.009302422 max resid 0.0738758
## Run 310 stress 0.1999554
## Run 311 stress 0.2006365
## Run 312 stress 0.1997217
## Run 313 stress 0.1991558
## Run 314 stress 0.1995132
## Run 315 stress 0.1995245
## Run 316 stress 0.19897
## Run 317 stress 0.2042452
## Run 318 stress 0.20049
## Run 319 stress 0.1989376
## Run 320 stress 0.1998948
## Run 321 stress 0.198354
## ... Procrustes: rmse 0.01368267 max resid 0.1316551
## Run 322 stress 0.2111917
## Run 323 stress 0.1990312
## Run 324 stress 0.1987924
## Run 325 stress 0.1989689
## Run 326 stress 0.1985801
## ... Procrustes: rmse 0.009976855 max resid 0.09460619
## Run 327 stress 0.1981951
## ... Procrustes: rmse 0.003923216 max resid 0.03818792
## Run 328 stress 0.1996472
## Run 329 stress 0.1991396
## Run 330 stress 0.1987015
## Run 331 stress 0.1996109
## Run 332 stress 0.1984842
## ... Procrustes: rmse 0.02425436 max resid 0.1580872
## Run 333 stress 0.2023409
## Run 334 stress 0.1995046
## Run 335 stress 0.1988666
## Run 336 stress 0.199047
## Run 337 stress 0.1998431
## Run 338 stress 0.1994146
## Run 339 stress 0.1988181
## Run 340 stress 0.1995038
## Run 341 stress 0.1988185
## Run 342 stress 0.2001735
## Run 343 stress 0.1985714
## ... Procrustes: rmse 0.02771676 max resid 0.175535
## Run 344 stress 0.200404
## Run 345 stress 0.2001418
## Run 346 stress 0.1990133
## Run 347 stress 0.2098598
## Run 348 stress 0.2011859
## Run 349 stress 0.1994488
## Run 350 stress 0.2001073
## Run 351 stress 0.1992629
## Run 352 stress 0.1994929
## Run 353 stress 0.1987395
## Run 354 stress 0.1991216
## Run 355 stress 0.198851
## Run 356 stress 0.2006935
## Run 357 stress 0.2006023
## Run 358 stress 0.2035624
## Run 359 stress 0.1990507
## Run 360 stress 0.1991665
## Run 361 stress 0.1987671
## Run 362 stress 0.1990774
## Run 363 stress 0.1987597
## Run 364 stress 0.199083
## Run 365 stress 0.1988346
## Run 366 stress 0.1989488
## Run 367 stress 0.199447
## Run 368 stress 0.1987979
## Run 369 stress 0.1989522
## Run 370 stress 0.1988724
## Run 371 stress 0.1995322
## Run 372 stress 0.1997618
## Run 373 stress 0.1995449
## Run 374 stress 0.2039276
## Run 375 stress 0.1988875
## Run 376 stress 0.1989139
## Run 377 stress 0.1988629
## Run 378 stress 0.1991641
## Run 379 stress 0.2014512
## Run 380 stress 0.1993177
## Run 381 stress 0.2007105
## Run 382 stress 0.1987594
## Run 383 stress 0.1986776
## ... Procrustes: rmse 0.01709909 max resid 0.1509284
## Run 384 stress 0.1987493
## Run 385 stress 0.2004005
## Run 386 stress 0.1987199
## Run 387 stress 0.2001559
## Run 388 stress 0.1988465
## Run 389 stress 0.1993712
## Run 390 stress 0.1992557
## Run 391 stress 0.2021603
## Run 392 stress 0.1994653
## Run 393 stress 0.2006692
## Run 394 stress 0.2017573
## Run 395 stress 0.2014722
## Run 396 stress 0.1996045
## Run 397 stress 0.199415
## Run 398 stress 0.198957
## Run 399 stress 0.1993117
## Run 400 stress 0.1987988
## Run 401 stress 0.1995233
## Run 402 stress 0.1986416
## ... Procrustes: rmse 0.02807089 max resid 0.1720642
## Run 403 stress 0.1992167
## Run 404 stress 0.1989215
## Run 405 stress 0.1997375
## Run 406 stress 0.2023143
## Run 407 stress 0.199448
## Run 408 stress 0.1994795
## Run 409 stress 0.2003865
## Run 410 stress 0.1989869
## Run 411 stress 0.1981909
## ... Procrustes: rmse 0.009716545 max resid 0.1195261
## Run 412 stress 0.2125122
## Run 413 stress 0.1999992
## Run 414 stress 0.1995522
## Run 415 stress 0.2000838
## Run 416 stress 0.199048
## Run 417 stress 0.2036027
## Run 418 stress 0.2003422
## Run 419 stress 0.19952
## Run 420 stress 0.1994282
## Run 421 stress 0.1984488
## ... Procrustes: rmse 0.02422448 max resid 0.1539348
## Run 422 stress 0.2005261
## Run 423 stress 0.2097327
## Run 424 stress 0.1987367
## Run 425 stress 0.1988187
## Run 426 stress 0.1991144
## Run 427 stress 0.2015914
## Run 428 stress 0.1991714
## Run 429 stress 0.1991732
## Run 430 stress 0.1998744
## Run 431 stress 0.1984901
## ... Procrustes: rmse 0.01917336 max resid 0.1580944
## Run 432 stress 0.2004334
## Run 433 stress 0.1992958
## Run 434 stress 0.1993535
## Run 435 stress 0.1991477
## Run 436 stress 0.1987315
## Run 437 stress 0.1987568
## Run 438 stress 0.1985556
## ... Procrustes: rmse 0.03136378 max resid 0.1810209
## Run 439 stress 0.1988589
## Run 440 stress 0.1989637
## Run 441 stress 0.2050953
## Run 442 stress 0.2027025
## Run 443 stress 0.2089542
## Run 444 stress 0.1999135
## Run 445 stress 0.2008746
## Run 446 stress 0.198991
## Run 447 stress 0.2001345
## Run 448 stress 0.1998776
## Run 449 stress 0.1988777
## Run 450 stress 0.1993712
## Run 451 stress 0.1991427
## Run 452 stress 0.2041956
## Run 453 stress 0.1994188
## Run 454 stress 0.1990918
## Run 455 stress 0.1995845
## Run 456 stress 0.1988127
## Run 457 stress 0.1988169
## Run 458 stress 0.1990798
## Run 459 stress 0.1982514
## ... Procrustes: rmse 0.009713179 max resid 0.08419412
## Run 460 stress 0.1988023
## Run 461 stress 0.2019274
## Run 462 stress 0.1982285
## ... Procrustes: rmse 0.007096129 max resid 0.0760459
## Run 463 stress 0.2022692
## Run 464 stress 0.1992604
## Run 465 stress 0.2016204
## Run 466 stress 0.1991247
## Run 467 stress 0.1996803
## Run 468 stress 0.1988167
## Run 469 stress 0.1986829
## Run 470 stress 0.1989017
## Run 471 stress 0.1990971
## Run 472 stress 0.1994561
## Run 473 stress 0.1990434
## Run 474 stress 0.2005555
## Run 475 stress 0.1987404
## Run 476 stress 0.1984279
## ... Procrustes: rmse 0.0172824 max resid 0.1562923
## Run 477 stress 0.1993171
## Run 478 stress 0.1987681
## Run 479 stress 0.1995565
## Run 480 stress 0.2000812
## Run 481 stress 0.2004669
## Run 482 stress 0.2001843
## Run 483 stress 0.2034905
## Run 484 stress 0.204827
## Run 485 stress 0.2019202
## Run 486 stress 0.2001194
## Run 487 stress 0.1988735
## Run 488 stress 0.1991146
## Run 489 stress 0.1981928
## ... Procrustes: rmse 0.008777207 max resid 0.1049153
## Run 490 stress 0.1995277
## Run 491 stress 0.1994786
## Run 492 stress 0.1990049
## Run 493 stress 0.1992867
## Run 494 stress 0.1991512
## Run 495 stress 0.1996907
## Run 496 stress 0.2129489
## Run 497 stress 0.1985096
## ... Procrustes: rmse 0.01895001 max resid 0.1582403
## Run 498 stress 0.198735
## Run 499 stress 0.1985077
## ... Procrustes: rmse 0.02839557 max resid 0.1750745
## Run 500 stress 0.1998199
## Run 501 stress 0.1995561
## Run 502 stress 0.1988606
## Run 503 stress 0.1987059
## Run 504 stress 0.1991007
## Run 505 stress 0.2004613
## Run 506 stress 0.1982131
## ... Procrustes: rmse 0.007463779 max resid 0.07480915
## Run 507 stress 0.2001035
## Run 508 stress 0.1984597
## ... Procrustes: rmse 0.02466766 max resid 0.1528222
## Run 509 stress 0.199124
## Run 510 stress 0.2013815
## Run 511 stress 0.1988257
## Run 512 stress 0.2104593
## Run 513 stress 0.1985613
## ... Procrustes: rmse 0.01505891 max resid 0.164745
## Run 514 stress 0.1988329
## Run 515 stress 0.1991055
## Run 516 stress 0.1991572
## Run 517 stress 0.2005056
## Run 518 stress 0.1994111
## Run 519 stress 0.1992869
## Run 520 stress 0.1988205
## Run 521 stress 0.1991169
## Run 522 stress 0.1991187
## Run 523 stress 0.1990306
## Run 524 stress 0.1983499
## ... Procrustes: rmse 0.01290028 max resid 0.1578367
## Run 525 stress 0.2019406
## Run 526 stress 0.1997294
## Run 527 stress 0.2021036
## Run 528 stress 0.2019837
## Run 529 stress 0.1988129
## Run 530 stress 0.1994107
## Run 531 stress 0.1994248
## Run 532 stress 0.1997897
## Run 533 stress 0.1994177
## Run 534 stress 0.1990885
## Run 535 stress 0.1988511
## Run 536 stress 0.1994809
## Run 537 stress 0.2000615
## Run 538 stress 0.1986578
## ... Procrustes: rmse 0.02133622 max resid 0.1429667
## Run 539 stress 0.1994873
## Run 540 stress 0.1992043
## Run 541 stress 0.1990678
## Run 542 stress 0.1991494
## Run 543 stress 0.1994112
## Run 544 stress 0.1991839
## Run 545 stress 0.2021392
## Run 546 stress 0.2103871
## Run 547 stress 0.1997408
## Run 548 stress 0.1985865
## ... Procrustes: rmse 0.01365693 max resid 0.1117205
## Run 549 stress 0.1994755
## Run 550 stress 0.1997297
## Run 551 stress 0.1992705
## Run 552 stress 0.1997218
## Run 553 stress 0.1995235
## Run 554 stress 0.2000416
## Run 555 stress 0.1991655
## Run 556 stress 0.1985973
## ... Procrustes: rmse 0.02432441 max resid 0.1658321
## Run 557 stress 0.1991257
## Run 558 stress 0.2029465
## Run 559 stress 0.1985716
## ... Procrustes: rmse 0.02698565 max resid 0.1550709
## Run 560 stress 0.1999401
## Run 561 stress 0.199196
## Run 562 stress 0.1988656
## Run 563 stress 0.1983411
## ... Procrustes: rmse 0.0109429 max resid 0.1067581
## Run 564 stress 0.1995066
## Run 565 stress 0.1991208
## Run 566 stress 0.1996409
## Run 567 stress 0.1990871
## Run 568 stress 0.1998312
## Run 569 stress 0.1983755
## ... Procrustes: rmse 0.0119995 max resid 0.1095442
## Run 570 stress 0.1990456
## Run 571 stress 0.1987428
## Run 572 stress 0.1992767
## Run 573 stress 0.1998159
## Run 574 stress 0.1986949
## Run 575 stress 0.1984499
## ... Procrustes: rmse 0.01313792 max resid 0.07894952
## Run 576 stress 0.1997635
## Run 577 stress 0.199782
## Run 578 stress 0.1988652
## Run 579 stress 0.1984942
## ... Procrustes: rmse 0.0165744 max resid 0.1287838
## Run 580 stress 0.2010908
## Run 581 stress 0.1991592
## Run 582 stress 0.1993268
## Run 583 stress 0.1988704
## Run 584 stress 0.1988905
## Run 585 stress 0.199469
## Run 586 stress 0.1992051
## Run 587 stress 0.1988066
## Run 588 stress 0.2031138
## Run 589 stress 0.199457
## Run 590 stress 0.1998772
## Run 591 stress 0.200046
## Run 592 stress 0.1986697
## ... Procrustes: rmse 0.02720056 max resid 0.1599086
## Run 593 stress 0.2022846
## Run 594 stress 0.1992039
## Run 595 stress 0.1992659
## Run 596 stress 0.1987785
## Run 597 stress 0.2006462
## Run 598 stress 0.200756
## Run 599 stress 0.1987982
## Run 600 stress 0.1996164
## Run 601 stress 0.1991959
## Run 602 stress 0.2033178
## Run 603 stress 0.2124951
## Run 604 stress 0.1990801
## Run 605 stress 0.198664
## ... Procrustes: rmse 0.03190207 max resid 0.1785979
## Run 606 stress 0.1995834
## Run 607 stress 0.2004999
## Run 608 stress 0.1993911
## Run 609 stress 0.198828
## Run 610 stress 0.1994561
## Run 611 stress 0.1990533
## Run 612 stress 0.1987401
## Run 613 stress 0.1989587
## Run 614 stress 0.2025542
## Run 615 stress 0.1988655
## Run 616 stress 0.1986641
## ... Procrustes: rmse 0.02852966 max resid 0.1729425
## Run 617 stress 0.1999527
## Run 618 stress 0.2092633
## Run 619 stress 0.2131853
## Run 620 stress 0.2027208
## Run 621 stress 0.2004011
## Run 622 stress 0.1981792
## ... Procrustes: rmse 0.009181807 max resid 0.1118831
## Run 623 stress 0.1992327
## Run 624 stress 0.1987568
## Run 625 stress 0.1986394
## ... Procrustes: rmse 0.02423259 max resid 0.1571553
## Run 626 stress 0.1989245
## Run 627 stress 0.2004719
## Run 628 stress 0.1987331
## Run 629 stress 0.2027782
## Run 630 stress 0.1995355
## Run 631 stress 0.2067446
## Run 632 stress 0.19936
## Run 633 stress 0.2117479
## Run 634 stress 0.1996739
## Run 635 stress 0.1982354
## ... Procrustes: rmse 0.01268791 max resid 0.1356392
## Run 636 stress 0.1990496
## Run 637 stress 0.1998852
## Run 638 stress 0.2005458
## Run 639 stress 0.2014989
## Run 640 stress 0.2022878
## Run 641 stress 0.2008355
## Run 642 stress 0.2007382
## Run 643 stress 0.2140726
## Run 644 stress 0.1999479
## Run 645 stress 0.1990361
## Run 646 stress 0.1992282
## Run 647 stress 0.1985568
## ... Procrustes: rmse 0.00968623 max resid 0.1104974
## Run 648 stress 0.1986071
## ... Procrustes: rmse 0.02921395 max resid 0.1740166
## Run 649 stress 0.1988914
## Run 650 stress 0.2012077
## Run 651 stress 0.1988652
## Run 652 stress 0.1985736
## ... Procrustes: rmse 0.02769806 max resid 0.1762039
## Run 653 stress 0.1989664
## Run 654 stress 0.201255
## Run 655 stress 0.201529
## Run 656 stress 0.1994893
## Run 657 stress 0.1987896
## Run 658 stress 0.2025157
## Run 659 stress 0.199411
## Run 660 stress 0.213469
## Run 661 stress 0.1985851
## ... Procrustes: rmse 0.02782007 max resid 0.1588786
## Run 662 stress 0.1988203
## Run 663 stress 0.1993199
## Run 664 stress 0.2031762
## Run 665 stress 0.2144611
## Run 666 stress 0.1988158
## Run 667 stress 0.1991153
## Run 668 stress 0.1994488
## Run 669 stress 0.1986807
## Run 670 stress 0.2105767
## Run 671 stress 0.2001535
## Run 672 stress 0.2050173
## Run 673 stress 0.2119223
## Run 674 stress 0.1985012
## ... Procrustes: rmse 0.01226937 max resid 0.1119371
## Run 675 stress 0.1981964
## ... Procrustes: rmse 0.004403934 max resid 0.04324853
## Run 676 stress 0.2011622
## Run 677 stress 0.1987966
## Run 678 stress 0.1990401
## Run 679 stress 0.1995558
## Run 680 stress 0.1990144
## Run 681 stress 0.1987698
## Run 682 stress 0.1986707
## ... Procrustes: rmse 0.02235539 max resid 0.1584492
## Run 683 stress 0.2007967
## Run 684 stress 0.1994541
## Run 685 stress 0.1999412
## Run 686 stress 0.2017024
## Run 687 stress 0.1988688
## Run 688 stress 0.1990354
## Run 689 stress 0.1991802
## Run 690 stress 0.1987022
## Run 691 stress 0.1982971
## ... Procrustes: rmse 0.006429099 max resid 0.07415036
## Run 692 stress 0.199029
## Run 693 stress 0.2001554
## Run 694 stress 0.1994725
## Run 695 stress 0.198808
## Run 696 stress 0.1992504
## Run 697 stress 0.1987644
## Run 698 stress 0.1986481
## ... Procrustes: rmse 0.02093763 max resid 0.2047165
## Run 699 stress 0.1995662
## Run 700 stress 0.1989153
## Run 701 stress 0.2005332
## Run 702 stress 0.1985695
## ... Procrustes: rmse 0.02526061 max resid 0.1530252
## Run 703 stress 0.198465
## ... Procrustes: rmse 0.02500176 max resid 0.1606646
## Run 704 stress 0.1987827
## Run 705 stress 0.1995575
## Run 706 stress 0.2007147
## Run 707 stress 0.1993846
## Run 708 stress 0.1988379
## Run 709 stress 0.1994774
## Run 710 stress 0.1988858
## Run 711 stress 0.1994234
## Run 712 stress 0.199241
## Run 713 stress 0.2010624
## Run 714 stress 0.1998403
## Run 715 stress 0.1993309
## Run 716 stress 0.1986153
## ... Procrustes: rmse 0.009215388 max resid 0.0754745
## Run 717 stress 0.1999681
## Run 718 stress 0.2019334
## Run 719 stress 0.1993625
## Run 720 stress 0.1991341
## Run 721 stress 0.1981933
## ... Procrustes: rmse 0.003896269 max resid 0.03968538
## Run 722 stress 0.1995463
## Run 723 stress 0.1990811
## Run 724 stress 0.1993889
## Run 725 stress 0.1989655
## Run 726 stress 0.1985904
## ... Procrustes: rmse 0.02486578 max resid 0.1682948
## Run 727 stress 0.1988484
## Run 728 stress 0.2006049
## Run 729 stress 0.1984839
## ... Procrustes: rmse 0.0114134 max resid 0.1109895
## Run 730 stress 0.2022474
## Run 731 stress 0.1992664
## Run 732 stress 0.1994694
## Run 733 stress 0.1988586
## Run 734 stress 0.2003405
## Run 735 stress 0.1992508
## Run 736 stress 0.1987902
## Run 737 stress 0.1986508
## ... Procrustes: rmse 0.02074596 max resid 0.148598
## Run 738 stress 0.2004612
## Run 739 stress 0.1987655
## Run 740 stress 0.1996411
## Run 741 stress 0.2014885
## Run 742 stress 0.1991255
## Run 743 stress 0.1994068
## Run 744 stress 0.1999323
## Run 745 stress 0.1981926
## ... Procrustes: rmse 0.004459676 max resid 0.04486127
## Run 746 stress 0.2002803
## Run 747 stress 0.2111272
## Run 748 stress 0.2032845
## Run 749 stress 0.1993874
## Run 750 stress 0.1988899
## Run 751 stress 0.200038
## Run 752 stress 0.1995965
## Run 753 stress 0.2115999
## Run 754 stress 0.1991943
## Run 755 stress 0.1984946
## ... Procrustes: rmse 0.0125227 max resid 0.1111627
## Run 756 stress 0.2014649
## Run 757 stress 0.198625
## ... Procrustes: rmse 0.02795817 max resid 0.1746954
## Run 758 stress 0.1994849
## Run 759 stress 0.1986523
## ... Procrustes: rmse 0.02211832 max resid 0.1575127
## Run 760 stress 0.1986573
## ... Procrustes: rmse 0.02824112 max resid 0.1729606
## Run 761 stress 0.1989098
## Run 762 stress 0.198562
## ... Procrustes: rmse 0.01029435 max resid 0.107689
## Run 763 stress 0.1991184
## Run 764 stress 0.1989355
## Run 765 stress 0.2112124
## Run 766 stress 0.2005146
## Run 767 stress 0.1989957
## Run 768 stress 0.198721
## Run 769 stress 0.1988732
## Run 770 stress 0.2130958
## Run 771 stress 0.1992434
## Run 772 stress 0.1988333
## Run 773 stress 0.1990951
## Run 774 stress 0.1988314
## Run 775 stress 0.1993185
## Run 776 stress 0.200409
## Run 777 stress 0.1989615
## Run 778 stress 0.2028477
## Run 779 stress 0.2032861
## Run 780 stress 0.199781
## Run 781 stress 0.199642
## Run 782 stress 0.1985033
## ... Procrustes: rmse 0.02845888 max resid 0.1762257
## Run 783 stress 0.1988356
## Run 784 stress 0.20062
## Run 785 stress 0.1987811
## Run 786 stress 0.2016737
## Run 787 stress 0.1994149
## Run 788 stress 0.1995355
## Run 789 stress 0.1987326
## Run 790 stress 0.1997353
## Run 791 stress 0.198959
## Run 792 stress 0.1984464
## ... Procrustes: rmse 0.02406715 max resid 0.1539157
## Run 793 stress 0.1999386
## Run 794 stress 0.1989379
## Run 795 stress 0.1990611
## Run 796 stress 0.1992599
## Run 797 stress 0.1986337
## ... Procrustes: rmse 0.01948131 max resid 0.1556282
## Run 798 stress 0.1994949
## Run 799 stress 0.199191
## Run 800 stress 0.1993729
## Run 801 stress 0.199018
## Run 802 stress 0.1994484
## Run 803 stress 0.1990973
## Run 804 stress 0.2108997
## Run 805 stress 0.1990685
## Run 806 stress 0.200234
## Run 807 stress 0.1991688
## Run 808 stress 0.1986407
## ... Procrustes: rmse 0.02799559 max resid 0.1733818
## Run 809 stress 0.1987568
## Run 810 stress 0.1996852
## Run 811 stress 0.1995678
## Run 812 stress 0.1998169
## Run 813 stress 0.198851
## Run 814 stress 0.1989062
## Run 815 stress 0.1994196
## Run 816 stress 0.198596
## ... Procrustes: rmse 0.02433934 max resid 0.1672624
## Run 817 stress 0.1994766
## Run 818 stress 0.200297
## Run 819 stress 0.19846
## ... Procrustes: rmse 0.02371486 max resid 0.1579017
## Run 820 stress 0.1986242
## ... Procrustes: rmse 0.02550899 max resid 0.1695047
## Run 821 stress 0.1990899
## Run 822 stress 0.1998445
## Run 823 stress 0.1989534
## Run 824 stress 0.1984141
## ... Procrustes: rmse 0.0174925 max resid 0.1591853
## Run 825 stress 0.1990156
## Run 826 stress 0.1987622
## Run 827 stress 0.1984301
## ... Procrustes: rmse 0.01706636 max resid 0.1565368
## Run 828 stress 0.1985582
## ... Procrustes: rmse 0.01167725 max resid 0.107879
## Run 829 stress 0.1994016
## Run 830 stress 0.1989631
## Run 831 stress 0.1988324
## Run 832 stress 0.2019106
## Run 833 stress 0.1991546
## Run 834 stress 0.2104972
## Run 835 stress 0.1986295
## ... Procrustes: rmse 0.009452878 max resid 0.07541767
## Run 836 stress 0.2011298
## Run 837 stress 0.1994193
## Run 838 stress 0.1994823
## Run 839 stress 0.2004076
## Run 840 stress 0.2005299
## Run 841 stress 0.2000567
## Run 842 stress 0.2003906
## Run 843 stress 0.1991622
## Run 844 stress 0.2021023
## Run 845 stress 0.1997408
## Run 846 stress 0.1986785
## ... Procrustes: rmse 0.02837167 max resid 0.1810907
## Run 847 stress 0.198588
## ... Procrustes: rmse 0.03108561 max resid 0.1746463
## Run 848 stress 0.200949
## Run 849 stress 0.1993652
## Run 850 stress 0.1987476
## Run 851 stress 0.1988207
## Run 852 stress 0.1997441
## Run 853 stress 0.198857
## Run 854 stress 0.2017857
## Run 855 stress 0.1995658
## Run 856 stress 0.2025617
## Run 857 stress 0.198791
## Run 858 stress 0.207986
## Run 859 stress 0.1993686
## Run 860 stress 0.213011
## Run 861 stress 0.2120092
## Run 862 stress 0.1994225
## Run 863 stress 0.2003732
## Run 864 stress 0.1987769
## Run 865 stress 0.1988656
## Run 866 stress 0.2001307
## Run 867 stress 0.1991365
## Run 868 stress 0.1991215
## Run 869 stress 0.1997306
## Run 870 stress 0.2028758
## Run 871 stress 0.1983106
## ... Procrustes: rmse 0.01309562 max resid 0.1125374
## Run 872 stress 0.1987979
## Run 873 stress 0.1990831
## Run 874 stress 0.1992495
## Run 875 stress 0.1997797
## Run 876 stress 0.1993077
## Run 877 stress 0.2001754
## Run 878 stress 0.1994319
## Run 879 stress 0.2005045
## Run 880 stress 0.1993236
## Run 881 stress 0.1989813
## Run 882 stress 0.201393
## Run 883 stress 0.1991007
## Run 884 stress 0.1990911
## Run 885 stress 0.1998224
## Run 886 stress 0.1987039
## Run 887 stress 0.1991685
## Run 888 stress 0.2000889
## Run 889 stress 0.1987639
## Run 890 stress 0.1989057
## Run 891 stress 0.1989074
## Run 892 stress 0.198823
## Run 893 stress 0.2004835
## Run 894 stress 0.198646
## ... Procrustes: rmse 0.02461157 max resid 0.1568601
## Run 895 stress 0.2004781
## Run 896 stress 0.2083318
## Run 897 stress 0.2020172
## Run 898 stress 0.1982624
## ... Procrustes: rmse 0.01230872 max resid 0.125031
## Run 899 stress 0.2041222
## Run 900 stress 0.2000858
## Run 901 stress 0.200085
## Run 902 stress 0.1988254
## Run 903 stress 0.2006446
## Run 904 stress 0.1990407
## Run 905 stress 0.1981843
## ... Procrustes: rmse 0.01009757 max resid 0.1275913
## Run 906 stress 0.1991798
## Run 907 stress 0.1988251
## Run 908 stress 0.2004603
## Run 909 stress 0.1994006
## Run 910 stress 0.1992692
## Run 911 stress 0.1988032
## Run 912 stress 0.1986353
## ... Procrustes: rmse 0.009735516 max resid 0.07571546
## Run 913 stress 0.1985715
## ... Procrustes: rmse 0.02756762 max resid 0.1748809
## Run 914 stress 0.1995393
## Run 915 stress 0.1994765
## Run 916 stress 0.2033272
## Run 917 stress 0.199608
## Run 918 stress 0.1991224
## Run 919 stress 0.1992652
## Run 920 stress 0.2005492
## Run 921 stress 0.2023237
## Run 922 stress 0.199657
## Run 923 stress 0.1988214
## Run 924 stress 0.2004266
## Run 925 stress 0.1991057
## Run 926 stress 0.2002619
## Run 927 stress 0.1986847
## Run 928 stress 0.1987496
## Run 929 stress 0.1987897
## Run 930 stress 0.1995559
## Run 931 stress 0.1986094
## ... Procrustes: rmse 0.02505669 max resid 0.168552
## Run 932 stress 0.1989428
## Run 933 stress 0.1985614
## ... Procrustes: rmse 0.01706529 max resid 0.1956583
## Run 934 stress 0.2002252
## Run 935 stress 0.1991314
## Run 936 stress 0.1989866
## Run 937 stress 0.1984874
## ... Procrustes: rmse 0.01558231 max resid 0.1112598
## Run 938 stress 0.1985011
## ... Procrustes: rmse 0.02843922 max resid 0.1757816
## Run 939 stress 0.1987761
## Run 940 stress 0.1986537
## ... Procrustes: rmse 0.01691212 max resid 0.1959829
## Run 941 stress 0.1993712
## Run 942 stress 0.1998255
## Run 943 stress 0.2022499
## Run 944 stress 0.1991754
## Run 945 stress 0.2101679
## Run 946 stress 0.1995377
## Run 947 stress 0.1987356
## Run 948 stress 0.2046635
## Run 949 stress 0.1996995
## Run 950 stress 0.1987171
## Run 951 stress 0.1997209
## Run 952 stress 0.1986367
## ... Procrustes: rmse 0.0200757 max resid 0.1456154
## Run 953 stress 0.1989731
## Run 954 stress 0.1989317
## Run 955 stress 0.1981693
## ... New best solution
## ... Procrustes: rmse 0.009022896 max resid 0.1127499
## Run 956 stress 0.1993496
## Run 957 stress 0.2000322
## Run 958 stress 0.2112859
## Run 959 stress 0.1986457
## ... Procrustes: rmse 0.02110193 max resid 0.1576134
## Run 960 stress 0.1984356
## ... Procrustes: rmse 0.01601182 max resid 0.1598652
## Run 961 stress 0.2026566
## Run 962 stress 0.1988674
## Run 963 stress 0.1991225
## Run 964 stress 0.2108145
## Run 965 stress 0.2016556
## Run 966 stress 0.1990825
## Run 967 stress 0.199382
## Run 968 stress 0.1991362
## Run 969 stress 0.1990412
## Run 970 stress 0.1995995
## Run 971 stress 0.1986997
## Run 972 stress 0.1989791
## Run 973 stress 0.1982418
## ... Procrustes: rmse 0.007015814 max resid 0.07418104
## Run 974 stress 0.2004303
## Run 975 stress 0.1995628
## Run 976 stress 0.1990294
## Run 977 stress 0.1992065
## Run 978 stress 0.1998431
## Run 979 stress 0.1994275
## Run 980 stress 0.2121839
## Run 981 stress 0.199664
## Run 982 stress 0.2011004
## Run 983 stress 0.1994721
## Run 984 stress 0.1996044
## Run 985 stress 0.1988742
## Run 986 stress 0.1994548
## Run 987 stress 0.1991597
## Run 988 stress 0.1981855
## ... Procrustes: rmse 0.002194565 max resid 0.02416401
## Run 989 stress 0.2122768
## Run 990 stress 0.2109429
## Run 991 stress 0.1995916
## Run 992 stress 0.1986419
## ... Procrustes: rmse 0.0183426 max resid 0.2158596
## Run 993 stress 0.1987934
## Run 994 stress 0.1997989
## Run 995 stress 0.2126021
## Run 996 stress 0.1988633
## Run 997 stress 0.1985747
## ... Procrustes: rmse 0.0291126 max resid 0.1741584
## Run 998 stress 0.198513
## ... Procrustes: rmse 0.01391815 max resid 0.1095238
## Run 999 stress 0.1994288
## Run 1000 stress 0.198791
## *** No convergence -- monoMDS stopping criteria:
## 127: no. of iterations >= maxit
## 873: stress ratio > sratmax
#prepare NMDS points and values for graphing
md.nmds.points <- md.nmds$points %>% # take the NMDS points
data.frame() %>% # make a data frame
rownames_to_column(var = "mysamples") %>%
mutate(SampleID = mysamples) %>%
column_to_rownames(var = "mysamples")
meta <- input_all$map_loaded %>%
rownames_to_column("SampleID")
md.nmds.metadata <- inner_join(x=md.nmds.points, y=meta, by = "SampleID") %>%
group_by(Life_Stage_Simplified) %>%
ungroup() %>%
dplyr::select(SampleID, MDS1, MDS2, everything()) # select these columns only
md.nmds.metadata.unq <- md.nmds.metadata %>%
dplyr::select(Life_Stage_Simplified, Type) %>%
unique()
md.nmds.metadata.unq$Life_Stage_Simplified <-
factor(md.nmds.metadata.unq$Life_Stage_Simplified,
levels=samptype_list)
# get stress value
stress.md = paste("stress =", round(md.nmds$stress, digits = 4))
# get NMDS hulls
group.chulls <- plyr::ddply(md.nmds.metadata, "Life_Stage_Simplified",
function(df) df[chull(df$MDS1, df$MDS2), ])
# Plot NMDS
nmds.plot <- ggplot(md.nmds.metadata, aes(x = MDS1, y = MDS2, shape = Type,
colour = Life_Stage_Simplified)) +
geom_polygon(data = group.chulls,
aes(fill=Life_Stage_Simplified),
alpha = 0.15, linetype=0) +
geom_point(size=3) +
annotate("text", x = Inf, y = Inf, label = stress.md, hjust = 1, vjust = 1) +
ggtitle("C") +
labs(colour="Sample Type") +
scale_colour_manual(values=samptype_cols) +
scale_fill_manual(values=samptype_cols) +
guides(fill=F, color = guide_legend(order = 1, override.aes=list(shape=15, size=6)),
shape = guide_legend(order = 0)) +
scale_shape_discrete(name = "", labels = c("Environment", "Toad")) +
theme(axis.title = element_text(size=16),
axis.text = element_text(size=14),
legend.title = element_text(size=16),
legend.text = element_text(size=14),
plot.title = element_text(size = 30))
# Make dispersion graph
# PERMDisp on sample type, strata=site
#create B-C species abundance matrix
input_all_rel <- convert_to_relative_abundances(input_all)
BC <- calc_dm(input_all_rel$data_loaded, method = "bray")
life.disper.lump <- betadisper(BC, input_all$map_loaded$Life_Stage_Simplified)
permutest(life.disper.lump)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 7 0.30517 0.043596 7.9241 999 0.001 ***
## Residuals 164 0.90228 0.005502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PermDisp_life_lump <- TukeyHSD(life.disper.lump) #Tukey post-hoc
kable(PermDisp_life_lump$group) %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| diff | lwr | upr | p adj | |
|---|---|---|---|---|
| Sediment-Soil | 0.1044081 | -0.0347484 | 0.2435645 | 0.2979445 |
| Water-Soil | 0.1570590 | 0.0160511 | 0.2980668 | 0.0176037 |
| Eggs-Soil | 0.1067868 | -0.0402240 | 0.2537976 | 0.3393009 |
| Tadpole-Soil | 0.1382782 | 0.0038326 | 0.2727239 | 0.0389765 |
| Metamorph-Soil | 0.0921875 | -0.0469690 | 0.2313439 | 0.4623355 |
| Subadult-Soil | 0.0043954 | -0.1619284 | 0.1707191 | 1.0000000 |
| Adult-Soil | 0.0281439 | -0.1151446 | 0.1714323 | 0.9988014 |
| Water-Sediment | 0.0526509 | -0.0156736 | 0.1209753 | 0.2652093 |
| Eggs-Sediment | 0.0023787 | -0.0776038 | 0.0823613 | 1.0000000 |
| Tadpole-Sediment | 0.0338702 | -0.0196150 | 0.0873554 | 0.5226472 |
| Metamorph-Sediment | -0.0122206 | -0.0766375 | 0.0521963 | 0.9990429 |
| Subadult-Sediment | -0.1000127 | -0.2115861 | 0.0115607 | 0.1147460 |
| Adult-Sediment | -0.0762642 | -0.1491792 | -0.0033492 | 0.0333328 |
| Eggs-Water | -0.0502721 | -0.1334340 | 0.0328897 | 0.5828757 |
| Tadpole-Water | -0.0187807 | -0.0769130 | 0.0393515 | 0.9750383 |
| Metamorph-Water | -0.0648715 | -0.1331959 | 0.0034530 | 0.0761105 |
| Subadult-Water | -0.1526636 | -0.2665377 | -0.0387895 | 0.0015396 |
| Adult-Water | -0.1289151 | -0.2053042 | -0.0525260 | 0.0000174 |
| Tadpole-Eggs | 0.0314914 | -0.0399812 | 0.1029640 | 0.8768006 |
| Metamorph-Eggs | -0.0145993 | -0.0945819 | 0.0653832 | 0.9992557 |
| Subadult-Eggs | -0.1023915 | -0.2236197 | 0.0188368 | 0.1658808 |
| Adult-Eggs | -0.0786430 | -0.1656157 | 0.0083298 | 0.1082157 |
| Metamorph-Tadpole | -0.0460908 | -0.0995760 | 0.0073945 | 0.1473621 |
| Subadult-Tadpole | -0.1338829 | -0.2395226 | -0.0282432 | 0.0035487 |
| Adult-Tadpole | -0.1101344 | -0.1735987 | -0.0466700 | 0.0000088 |
| Subadult-Metamorph | -0.0877921 | -0.1993655 | 0.0237812 | 0.2406801 |
| Adult-Metamorph | -0.0640436 | -0.1369586 | 0.0088714 | 0.1309603 |
| Adult-Subadult | 0.0237485 | -0.0929377 | 0.1404347 | 0.9984928 |
# write.table(PermDisp_life_lump$group,
# paste0(work_dir,"Figs_Tables/PermDisp_life_lump.txt"), sep="\t")
Disper1 <- data.frame(life.disper.lump$distances)
colnames(Disper1) <- "dispers"
Disper2 <- Disper1 %>%
rownames_to_column("SampleID") %>%
inner_join(y=meta,
by="SampleID")
disper_bars <- ggplot(data=Disper2, aes(x=Life_Stage_Simplified, y=dispers)) +
geom_boxplot(aes(fill=Life_Stage_Simplified), alpha=0.6,outlier.shape=NA) +
ylab("Dispersion") + xlab("") + ggtitle("B") +
geom_point(aes(fill=Life_Stage_Simplified), size=3, shape=21, position=position_jitter(
width=0.25
)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
axis.title = element_text(size=16),
axis.text = element_text(size=14),
legend.title = element_text(size=16),
legend.text = element_text(size=14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(size = 30)) +
labs(fill="Sample Type") +
scale_fill_manual(values=samptype_cols)
legend <- arrangeGrob(get_legend(disper_bars))
disper_bars_noleg <- disper_bars + theme(legend.position = "none")
# Put all into paneled graph with annotations and shared legend
diversity_fig <- ggpubr::ggarrange(arrangeGrob(OTUrich, disper_bars_noleg),
legend, nmds.plot, ncol=3,
widths=c(4.5,1.5,9))
diversity_fig
ggsave(paste0(work_dir,"Figs_Tables/FungDivFig.png"), diversity_fig, dpi=800,
width=16, height=7)Permanova for the above plot
# PERMANOVA on life stage, strata=site
Perm <- adonis(BC ~ input_all$map_loaded$Life_Stage_Simplified,
strata = input_all$map_loaded$Location)
kable(Perm$aov.tab) %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| input_all\(map_loaded\)Life_Stage_Simplified | 7 | 8.634371 | 1.233482 | 3.118215 | 0.1174611 | 0.001 |
| Residuals | 164 | 64.873969 | 0.395573 | NA | 0.8825389 | NA |
| Total | 171 | 73.508340 | NA | NA | 1.0000000 | NA |
# also show the pairwise PERMANOVA
pw_perm <- read.delim(paste0(work_dir,"Figs_Tables/PW_PERM.txt"))
kable(pw_perm) %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| X1 | X2 | R2 | pval | pvalBon | pvalFDR |
|---|---|---|---|---|---|
| Water | Sediment | 0.0971807 | 0.001 | 0.091 | 0.005 |
| Water | Eggs.11.15 | 0.0811505 | 0.001 | 0.091 | 0.005 |
| Water | Metamorph.40.46 | 0.0489309 | 0.001 | 0.091 | 0.005 |
| Water | Subadult | 0.0833049 | 0.001 | 0.091 | 0.005 |
| Water | Adult | 0.0987144 | 0.001 | 0.091 | 0.005 |
| Sediment | Soil | 0.0899016 | 0.001 | 0.091 | 0.005 |
| Sediment | Eggs.11.15 | 0.0821637 | 0.001 | 0.091 | 0.005 |
| Sediment | Tadpole.25.27 | 0.0999341 | 0.001 | 0.091 | 0.005 |
| Sediment | Tadpole.27.29 | 0.1020966 | 0.001 | 0.091 | 0.005 |
| Sediment | Tadpole.29.31 | 0.0919525 | 0.001 | 0.091 | 0.005 |
| Sediment | Tadpole.31.35 | 0.0855843 | 0.001 | 0.091 | 0.005 |
| Sediment | Metamorph.40.46 | 0.1337229 | 0.001 | 0.091 | 0.005 |
| Sediment | Subadult | 0.1161327 | 0.001 | 0.091 | 0.005 |
| Sediment | Adult | 0.1504275 | 0.001 | 0.091 | 0.005 |
| Soil | Tadpole.27.29 | 0.0702598 | 0.001 | 0.091 | 0.005 |
| Soil | Metamorph.40.46 | 0.0871272 | 0.001 | 0.091 | 0.005 |
| Soil | Adult | 0.1516369 | 0.001 | 0.091 | 0.005 |
| Eggs.11.15 | Tadpole.25.27 | 0.0884647 | 0.001 | 0.091 | 0.005 |
| Eggs.11.15 | Tadpole.27.29 | 0.0750862 | 0.001 | 0.091 | 0.005 |
| Eggs.11.15 | Metamorph.40.46 | 0.0959784 | 0.001 | 0.091 | 0.005 |
| Eggs.11.15 | Subadult | 0.1396777 | 0.001 | 0.091 | 0.005 |
| Eggs.11.15 | Adult | 0.1074929 | 0.001 | 0.091 | 0.005 |
| Tadpole.25.27 | Adult | 0.1082529 | 0.001 | 0.091 | 0.005 |
| Tadpole.27.29 | Metamorph.40.46 | 0.0405703 | 0.001 | 0.091 | 0.005 |
| Tadpole.27.29 | Subadult | 0.0688705 | 0.001 | 0.091 | 0.005 |
| Tadpole.27.29 | Adult | 0.0907215 | 0.001 | 0.091 | 0.005 |
| Tadpole.29.31 | Subadult | 0.1499271 | 0.001 | 0.091 | 0.005 |
| Tadpole.29.31 | Adult | 0.1046980 | 0.001 | 0.091 | 0.005 |
| Tadpole.31.35 | Adult | 0.0929699 | 0.001 | 0.091 | 0.005 |
| Tadpole.36.39 | Adult | 0.1088128 | 0.001 | 0.091 | 0.005 |
| Metamorph.40.46 | Subadult | 0.0825493 | 0.001 | 0.091 | 0.005 |
| Metamorph.40.46 | Adult | 0.1116956 | 0.001 | 0.091 | 0.005 |
| Subadult | Adult | 0.1024093 | 0.001 | 0.091 | 0.005 |
Q2: Are there transient versus symbiotic OTU’s on boreal toad skin?
Which fungi are more likely to be symbionts in that they are associated with the host as opposed to transiently moving through the host system?
To answer this question, I’m using ANCOM, an indicator species analysis for microbial sequence datasets. - Uses log-ratio transformation to somewhat account for compositionality of the data, which helps decrease the false discovery rate. - Uses K-W and wilcoxon tests to determine significance. - Seems to be more stringent than other methods, like LEfSe or macroecology indicator species analysis, in terms of how many OTU’s pass the thresholds. - Does not need rarefied data, but I am using my rarefied data here. - ANCOM is running a bunch of pairwise tests and the W stat is how many of those pairwise tests passed the “significance” threshold set in the initial script when I ran it.
Run ANCOM
# Ancom needs Sample.ID to be the header starting at A1, which means changing it in both files and transposing the taxa table.
# export_taxa_table(input = input_all, "taxa_table_ancom.txt", "map_ancom.txt")
# Did changes in excel, and changed OTU ID's to have the taxonomy and OTU # (for future reference)
# re-read in clean files
OTU_ancom <- read.delim(paste0(work_dir,"02_hypothesisbasedanalsysis/otu_ancom.txt"),
sep = "\t", header = T)
map_ancom <- read.delim(paste0(work_dir,"02_hypothesisbasedanalsysis/map_ancom.txt"),
sep = "\t", header = T)
# source the code, downloaded from: https://github.com/mech3132/MSC_code/blob/master/ANCOM_updated_code_MYC.R
source(paste0(work_dir,"02_hypothesisbasedanalsysis/ANCOM_updated_code_MYC.R"))
# 1) using just sample type as life or env
# run ancom with parameters for no repeated data, no adjustment (so it'll use a Wilcoxon test)
comparison_test_lifeenv <- ANCOM.main.myc(OTUdat = OTU_ancom,
Vardat = map_ancom,
adjusted = F,
repeated = F,
main.var = "Type",
adj.formula = NULL,
repeat.var = NULL,
longitudinal = F,
random.formula = NULL,
multcorr = 2,
sig = 0.05,
prev.cut = 0.90)
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Calculating critical statistic values..."
## [1] "Generating plots..."
## [1] "Done"
# This tells you if your data has a lot of zeroes/skew, as most microbial data do.
# everything below the red line gets removed.
comparison_test_lifeenv$PLot.zeroes
# This tells you which (colored) OTU's have a high effect size and significance
comparison_test_lifeenv$Plot.volcano
comparison_test_lifeenv$Plot.volcano.0.9
# this helps you determine if the critical w was chosen correctly. We want the valley between a bimodal distribution, in theory.
comparison_test_lifeenv$Plot.critical.w 
# I think W=~20 actually would be better, especially with the clumping of the points in the previous plots. I'll take this into account in the write out file.
# write.table(comparison_test_lifeenv$W.taxa, paste0(work_dir,"02_hypothesisbasedanalsysis/indic_taxabunds_lifeenv.txt"), sep="\t", quote=F)
# 2) ANCOM so that each life stage gets compared individually to each environmental sample type. Caveat: this uses a Kruskall-Wallis test so it's not a pairwise comparison (i.e. it won't identify OTU's that are different between e.g. water and tadpole)
comparison_test_lifestages <- ANCOM.main.myc(OTUdat = OTU_ancom,
Vardat = map_ancom,
adjusted = F,
repeated = F,
main.var = "Life_Stage_Simplified",
adj.formula = NULL,
repeat.var = NULL,
longitudinal = F,
random.formula = NULL,
multcorr = 2,
sig = 0.05,
prev.cut = 0.90)
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Calculating critical statistic values..."
## [1] "Generating plots..."
## [1] "Done"
# write.table(comparison_test_lifestages$W.taxa, paste0(work_dir,"02_hypothesisbasedanalsysis/indic_taxabundsnew.txt"), sep="\t", quote=F)Graph Volcano plot
#make a volcano plot showing the associated microbes
otu.names <- comparison_test_lifeenv$W.taxa$otu.names
W_stat <- comparison_test_lifeenv$W.taxa$W_stat
W_f <- comparison_test_lifeenv$W.taxa$W_f
W_frame <- data.frame(otu.names,W_stat,W_f,row.names=NULL)
# in excel, add column for whether they are transient, fac, obl; read in that file
# write.table(W_frame, paste0(work_dir,"Figs_Tables/Volcano_plot/volcano_plot.txt"), sep="\t")
W_frame_V2 <- read.delim(paste0(work_dir,"Figs_Tables/Volcano_plot/volcano_plot_V2.txt"))
26-(26-18)/2 # want to put interept between 18 and 26 because everything 26+ is significant and 18 and less is not
## [1] 22
# figure out order of ecol for legend
Ecol_list <- c("Obligate", "Facultative", "Transient", "Transient, not significant")
W_frame_V2$Ecol <- factor(W_frame_V2$Ecol, levels=Ecol_list)
# graph volcano plot
volcano.plot <- ggplot(W_frame_V2) +
geom_point(aes(x=W_f, y=W_stat, color=Ecol)) +
xlab("CLR F statistic (Effect size)") +
ylab("W statistic") +
labs(color = "") +
scale_color_manual(values=c("#0072B2","#56B4E9", "#CC79A7", "#000000"))+
geom_hline(mapping=aes(yintercept=22), col="red", lty=2)
volcano.plot
# ggsave(paste0(work_dir,"Figs_Tables/Volcano_plot/volcano_plot_color.png"),
# volcano.plot, dpi=400, width=7, height=5)Graph relative abundances for significant fungal OTU’s
I used excel to clean up the output file “indic_taxabunds_lifeenv.txt” (ANCOM output file) into “indic_taxabunds.txt” which has the information I need moving forward.
# read in cleaned up file with relative abundances
indic_abunds <- read.delim(paste0(work_dir,"02_hypothesisbasedanalsysis/indic_taxabunds.txt"), header=T) %>%
gather(key="Sample.Type", value="Taxon.Relative.Abundance", -c("OTU","Ecol")) %>%
mutate(Type=case_when(Sample.Type=="Soil" ~ "Environment",
Sample.Type=="Sediment" ~ "Environment",
Sample.Type=="Water" ~ "Environment",
Sample.Type=="Eggs" ~ "Toad",
Sample.Type=="Tadpole" ~ "Toad",
Sample.Type=="Metamorph" ~ "Toad",
Sample.Type=="Subadult" ~ "Toad",
Sample.Type=="Adult" ~ "Toad"))
indic_abunds$Sample.Type <-
factor(indic_abunds$Sample.Type,
levels=c("Soil", "Sediment", "Water", "Eggs", "Tadpole",
"Metamorph", "Subadult", "Adult"))
# make three sets of faceted graphs and put on a plot together
obl_abunds <- dplyr::filter(indic_abunds, Ecol == "Obligate")
fac_abunds <- dplyr::filter(indic_abunds, Ecol=="Facultative")
tran_abunds <- dplyr::filter(indic_abunds, Ecol=="Transient")
Obl_g <- ggplot(obl_abunds, aes(x=Sample.Type, y=Taxon.Relative.Abundance,
group=OTU, fill=Type)) +
geom_bar(stat="identity") +
xlab("") + ylab("") +
scale_fill_manual(values=c("#548235", "#C55A11"), name="Sample Type") +
theme(axis.text.x=element_text(angle=90, hjust=1, size=20),
axis.title = element_text(size=30),
axis.text.y=element_text(size=20),
strip.text = element_text(size = 20),
legend.position = "none") +
facet_wrap(~OTU, scales="free_y", ncol=3, labeller = label_wrap_gen(width=15))
Fac_g <- ggplot(fac_abunds, aes(x=Sample.Type, y=Taxon.Relative.Abundance,
group=OTU, fill=Type)) +
geom_bar(stat="identity") +
xlab("") + ylab("Taxon Relative Abundance") +
scale_fill_manual(values=c("#548235", "#C55A11"), name="Sample Type") +
theme(axis.text.x=element_text(angle=90, hjust=1, size=20),
axis.title = element_text(size=30),
axis.text.y=element_text(size=20),
strip.text = element_text(size = 20),
legend.title = element_text(size=30),
legend.text = element_text(size=20)) +
facet_wrap(~OTU, scales="free_y", ncol=3, labeller = label_wrap_gen(width=15))
Tran_g <- ggplot(tran_abunds, aes(x=Sample.Type, y=Taxon.Relative.Abundance,
group=OTU, fill=Type)) +
geom_bar(stat="identity") +
xlab("") + ylab("") +
scale_fill_manual(values=c("#548235", "#C55A11"), name="Sample Type") +
theme(axis.text.x=element_text(angle=90, hjust=1, size=20),
axis.title = element_text(size=30),
axis.text.y=element_text(size=20),
strip.text = element_text(size = 20),
legend.position = "none") +
facet_wrap(~OTU, scales="free_y", ncol=3, labeller = label_wrap_gen(width=15))
Obl_g; Fac_g; Tran_g
# pull out the legend (same one for all of the graphs)
legend_abund <- arrangeGrob(get_legend(Fac_g))
# ggsave(paste0(work_dir,"Figs_Tables/Symbionts_fig/legend_abund.png"), legend_abund, dpi=400,
# width=3, height=3)
Facg_noleg <- Fac_g + theme(legend.position = "none")
#save all these graphs and rearrange how I want them in powerpoint
# ggsave(paste0(work_dir,"Figs_Tables/Symbionts_fig/obl_graph.png"),
# Obl_g, dpi=400, width=12, height=4.5)
# ggsave(paste0(work_dir,"Figs_Tables/Symbionts_fig/fac_graph.png"),
# Facg_noleg, dpi=400, width=12, height=7)
# ggsave(paste0(work_dir,"Figs_Tables/Symbionts_fig/tran_graph.png"),
# Tran_g, dpi=400, width=12, height=4.5)
# Assembled these files into the published graph using powerpointMiscellaneous
Here, I’m making a graph showing where certain sample types fell out in the various processing steps.
drops <- read.delim(paste0(work_dir,"Figs_Tables/Metadata_dropouts_R.txt"))
#reorder the axes
all_sampletypes <- c("Upland_Soil", "Sediment", "Water", "Egg", "Tadpole", "Metamorph", "Subadult", "Adult", "Mouth", "Tadpole_Dead", "Glove", "Sterile_water", "NoTemplateControl")
drops$Life_Stage_Simplified <- factor(drops$Life_Stage_Simplified,
levels=all_sampletypes)
madetosteplist <- c("noDNA", "filtSeq", "filtR", "Final")
drops$Made_to_step <- factor(drops$Made_to_step, levels=madetosteplist)
#make graph
drop_graph <- ggplot(drops, aes(fill=Life_Stage_Simplified,
x=Made_to_step)) +
geom_bar(position="dodge", stat="count") +
labs(fill="Sample Type", x="Reasons for samples dropping out",
y="Count") +
scale_x_discrete(labels=c("PCR did not amplify target gene",
"Filtered during sequence processing
(low reads or quality)",
"Filtered during R analysis
(some sample types not useful
for answering hypotheses)",
"Final dataset")) +
theme(axis.text.x = element_text(angle=0, hjust=0.5, size=14),
axis.title = element_text(size=16),
axis.text.y = element_text(size=14),
legend.text = element_text(size=14),
legend.title = element_text(size=16))
drop_graph
# ggsave("drop_graph.jpg", drop_graph, width=21, height=10)